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Ion  thruster  plume  interaction  has  been  studied  extensively  in  recent  years.  While 
most  existing  plume  models  have  focused  on  charge-exchange  ion  interactions  with  the 
spacecraft  and/or  plume-induced  contamination,  few  studies  have  examined  the  detailed 
physics  near  the  thruster  exit.  In  particular,  the  ion  beam  neutralization  process  and  the 
characteristics  of  the  neutralizing  electrons  are  not  well  understood.  This  paper  presents 
a  full-PIC  model  for  the  near-thruster  plume  of  a  single  thruster.  A  multi-domain  model, 
which  splits  the  domain  into  a  near-cathode  region  and  a  near-plume  region,  is  used  to 
obtain  velocity  distribution  of  the  neutralizing  electrons.  This  simulation  is  compared  with 
one  that  assumes  a  pre-neutralized  beam  and  another  that  uses  a  floating  cathode  potential. 


I.  Introduction 

Numerical  modeling  has  been  used  extensively  to  study  the  interaction  of  electric  thruster  plume  with 
spacecraft.  Most  codes  simplify  the  plasma  dynamics  by  tracking  only  the  heavy  particles  (ions  and  neu¬ 
trals),  while  making  assumptions  about  the  distribution  of  the  electrons.  The  mass  of  the  electrons,  and  thus 
the  difference  in  time  scales  of  the  various  species,  makes  a  full  Particle-In-Cell  (PIC)1  algorithm  computa¬ 
tionally  challenging.  A  common  approach  is  to  assume  that  the  electron  distribution  follows  the  Boltzmann 
relationship,  based  on  user  supplied  values  of  reference  parameters.  A  typical  hybrid-PIC  algorithm  thus 
cannot  correctly  resolve  the  physics  in  the  near-thruster  region,  which  is  dominated  by  a  non-neutral  plasma 
and  an  interaction  of  neutralizing  electrons  with  the  beam.  No  simulation  models  are  currently  available  to 
investigate  the  near-thruster  plume,  and  the  ion  beam  neutralization  process  is  still  not  well  understood. 

Understanding  of  the  neutralization  process  will  become  even  more  important  for  electric  thruster  clusters 
that  are  being  considered  for  future  space  applications.  Such  a  cluster  system  may  use  a  single  cathode  to 
neutralize  ion  beams  emitted  from  multiple  thrusters.  Already,  several  cluster  configurations  were  tested 
experimentally  by  Beal  and  Hargus.2,3  However,  experimental  measurements  can  provide  only  a  limited 
amount  of  information  on  the  motion  of  the  electrons.  The  effectiveness  of  beam  neutralization  by  a  shared 
neutralizer  is  still  not  clear,  and  there  is  no  generally  accepted  optimal  configuration  for  electric  thruster 
clusters. 

This  paper  presents  results  from  numerical  modeling  of  neutralization  of  a  single  ion  thruster.  An  attempt 
has  been  made  to  improve  results  obtained  in  a  previous  work.4  Primarily,  the  effect  of  the  simplified 
cathode  model  on  electron  dynamics  is  investigated.  In  the  present  study,  a  multi-domain  formulation  was 
used  to  obtain  initial  distribution  of  electrons  near  the  cathode.  The  electrons  were  then  sampled  from  the 
distribution  list  and  then  introduced  into  the  ion  bema  domain.  This  approach  was  necessary,  because  the 
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main  simulation  mesh  was  not  fine  enough  to  resolve  the  plasma  parameters  in  the  high-density  region  near 
the  cathode  tip. 


II.  DRACO  ES-PIC  Code 

A  3D  plasma  simulation  code,  called  DRACO,  was  used  in  this  work.  DRACO  was  developed  within  the 
AFRL  COLISEUM  framework,  which  is  a  collection  of  modules  capable  of  modeling  the  dynamics  of  electric 
thruster  plumes  and  their  interactions  with  spacecraft  surfaces.5 

The  DRACO  module  tracks  particles  on  a  Cartesian  mesh,  which  has  been  overlaid  with  a  secondary 
tetrahedral  mesh.6  This  secondary  mesh  allows  DRACO  to  resolve  surface  geometries  with  detail  beyond 
the  standard  ” stair-case”  representation  attainable  on  Cartesian  grids.  Surface  definition  is  specified  using 
planar  cuts  of  interface  tetrahedrons. 

The  main  COLISEUM  package  contains  support  for  loading  of  triangular  meshes  from  input  files  using 
standard  formats  such  as  Ansys  or  Abaqus.  The  interface  mesh  is  generated  automatically  by  DRACO’s 
helper  module  called  VOLCAR.  The  actual  intersection  process  is  described  in  a  greater  detail  in  Ref. 7. 

The  interface  cuts  are  used  to  perform  particle  surface  interactions,  and,  depending  on  the  chosen  solver, 
to  obtain  the  plasma  potential,  <j).  The  plasma  potential  is  computed  from  the  Poisson’s  equation, 

V20  =  -^  (1) 

£o 

using  the  DADI8  scheme.  In  the  above  equation,  p  is  the  charge  density  of  the  particles,  C/m3,  and  £o 
is  the  permittivity  of  free  space,  8.854  x  10-12  F/m.  The  charge  density  is  computed  from  the  individual 
contributions  of  the  ions  and  the  electrons,  p  =  q(ni  —  ne),  where  is  the  number  density  of  the  ions  or 
electrons.  In  the  Particle-In-Cell  (PIC)  method,  the  number  density  is  obtained  by  coupling  the  particles 
with  the  grid  through  particle  shape  factors , 


nk=^2  wis  (xi  -  xk)  (2) 

i 

where  Xk  is  the  position  of  a  grid  node,  and  is  the  specific  weight  of  the  macroparticle.  In  this  work,  the 
shape  and  size  of  the  particles  was  identical  to  the  Cartesian  cell.  This  first-order  representation  reduces 
the  simulation  noise  associated  with  the  zeroeth-order  (point  particle)  model,  while  still  allowing  a  simple 
particle-mesh  weighing  algorithm.  The  electric  field,  E,  is  then  computed  from 

V0  =  —E  (3) 

using  the  standard  centered  finite-difference  method.  Particle  velocity  is  adjusted  according  to  the  Lorentz 
force, 

=  F  =  q(jtJ  +  vxB>j  (4) 


where  m 

=  particle  mass,  kg 

V 

=  particle  velocity,  m /: 

Q 

=  particle  charge,  C 

E 

=  electric  field,  V/m 

B 

=  magnetic  field,  T 

The  electro-static  (ES)  formulation,  implemented  by  DRACO,  assumes  that  dB/dt  =  0.  No  static  back¬ 
ground  field  was  used  in  the  current  simulation,  and  hence  the  force  acting  on  the  particles  was  simply 

dv  ->  -> 

m—  =  F  =  qE  (5) 


The  equation  of  motion  for  the  particles  is 


(6) 


This  equation  is  integrated  numerically  along  with  eq.  5  using  the  leapfrog  method  with  a  finite  timestep  At. 
Final  position  of  the  particles  is  checked  for  surface  interactions.  Particles  leaving  the  simulation  domain  are 
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either  removed  from  the  simulation,  or  are  reintroduced  according  to  prescribed  boundary  conditions.  New 
particles  are  introduced  by  sampling  particle  sources.  This  process  repeats  until  a  user  specified  condition, 
such  as  reaching  a  steady-state,  or  exceeding  a  maximum  number  of  timesteps,  is  satisfied. 

III.  Neutralization  Modeling 


A.  Thruster  Model 

The  ion  thruster  used  in  this  study  is  based  on  the  40cm  NASA  Evolutionary  Xenon  Thruster  (NEXT). 
Surface  definition  of  the  thruster  is  shown  in  figure  1.  This  mesh  was  generated  with  MSC.Patran.  It  should 
be  noted  that  a  dimensional  drawing  of  the  thruster  was  not  available  to  the  authors,  and  hence  the  thruster 
geometry  was  generated  by  collecting  data  from  several  relevant  sources.9,10 


(a)  x  —  y  plane 


x 


(b)  y  —  z  plane 


Figure  1.  Surface  mesh  of  the  ion  thruster.  Yellow  region  indicates  source  triangles  emitting  ions.  Electrons 
are  injected  from  the  small  blue  region  at  the  tip  of  the  cathode.  The  physical  curvature  of  the  ion  optics  was 
used  to  introduce  curvature  to  the  ion  beam.  The  normal  vectors  of  source  elements,  shown  by  red  arrows, 
roughly  indicate  the  initial  divergence  of  the  beam. 


COLISEUM  provides  a  general  support  for  surface  sources.  Various  particle  production  models  can  be 
associated  with  a  collection  of  surface  triangles.  Generally,  the  particles  are  introduced  with  a  mean  velocity 
in  the  direction  of  the  surface  normal  vector.  Physical  curvature  of  the  surface  will  result  in  a  divergence 
of  the  ion  beam.  This  feature  was  used  in  the  present  work,  as  can  be  seen  in  figure  1(b),  which  shows 
the  normal  vectors  of  the  surface  elements.  Curvature  of  the  surface  mesh  resulted  in  approximately  a  15° 
beam  divergence.  Beam  flatness  (ratio  of  current  density  between  the  centerline  and  the  edge)  was  adjusted 
by  biasing  the  mass  production  rate  of  the  source  elements,  according  to  current  density  measurements  of 
Soulas.11  The  thruster  operating  condition  was  set  to  1.2 A  of  beam  current,  with  ions  injected  with  an 
average  velocity  of  34,400m/s  (corresponding  to  3510s  ISP)  and  a  temperature  of  O.leV.12 

Several  assumptions  were  made  about  the  plume  dynamics.  First,  collisions  were  ignored.  The  mean  free 
path,  A m,  for  an  electron- ion  collision  can  be  estimated  from  Ref.13 


A 


m 


1 

na 


167T£Qm2V4 

neA 


(7) 


For  plasma  density,  n  ^  1015  m-3  and  electron  velocities,  v  ~  106  m/s,  Am  is  0(l)m.  This  length  is  similar  to 
the  characteristic  dimension  of  the  domain.  Coulomb  collisions  thus  do  not  play  a  significant  role.  Neutral- 
neutral  and  neutral- ion  collisions  are  expected  to  be  more  frequent,  leading  to  a  particle  scatter  and  creation 
of  charge- exchange  ions.  While  the  effect  of  these  collisions  is  significant  in  the  study  of  interaction  of  the 
plume  with  the  spacecraft,  this  work  concentrated  on  the  dynamics  of  the  electrons  and  their  interaction 
with  the  beam  ions. 

Second,  the  plume  was  assumed  to  consist  only  of  singly-charged  ions  and  electrons.  While  the  neutral 
density  near  the  thruster  exit  can  exceed  the  ion  density,  the  neutrals  interact  with  ions  exclusively  through 
collisions.  In  the  absence  of  collisions,  tracking  of  neutrals  would  only  slow  down  the  simulation.  Doubly- 
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charged  ions  were  not  included  since  their  actual  distribution  was  not  known.  Furthermore,  production  of 
doubly- charged  ions  is  not  desired  because  it  leads  to  faster  thruster  erosion. 

Finally,  the  thruster  was  assumed  to  be  a  perfect  conductor.  Any  absorbed  electrons  were  re-injected 
from  their  original  source  at  the  next  timestep.  Collection  of  an  electron  current  in  a  space  environment 
would  lead  to  a  decrease  in  the  thruster  potential,  thus  preventing  an  excessive  build-up  of  negative  charge. 
However,  the  current  version  of  DRACO  does  not  contain  functionality  to  model  this  behavior. 

B.  Dimensional  Scaling  and  Boundary  Conditions 

Interaction  of  the  electrons  with  the  ions,  and  their  containment  in  the  beam  was  examined  by  modeling  a 
pre-neutralized  beam.  This  beam  pre-neutralization  was  accomplished  by  injecting  both  ions  and  electrons 
from  the  optics.  Due  to  symmetry,  only  a  quarter  domain  was  simulated.  The  domain  is  shown  in  figure 
2(a).  The  cell  size  was  set  to  2  x  10-4  ^  A £>,  and  the  simulation  contained  50  x  50  x  90  cells. 


(a)  simulation  domain  (b)  charge  density 


Figure  2.  Simulation  domain  used  to  study  electron  dynamics  in  a  pre-neutralized  beam.  Right  plot  (b)  shows 
charge  density  at  steady-state  if  a  uniform  cell  size  of  2cm  was  used.  Large  cell  size  leads  to  formation  of  a 
virtual  anode. 


The  simulation  was  performed  on  a  thruster  scaled-down  by  a  factor  of  100:1.  Plasma  density  at  the 
thruster  exit  was  retained  by  decreasing  the  emission  current  by  10,000  (100  x  100).  This  scaling  was 
necessary,  since  resolving  the  Debye  length  on  the  full-sized  domain  would  require  a  numerically  excessive 
number  of  computational  nodes  (over  1  billion).  Total  number  of  nodes  could  be  reduced  by  increasing  the 
simulation  cell  size  to  about  100A^.  However,  as  figure  2(b)  shows,  doing  so  results  in  a  development  of 
a  virtual  anode14  at  the  thruster  exit,  despite  the  thruster  injecting  equal  electron  and  ion  currents.  The 
virtual  anode  develops  since  the  PIC  formulation  replaces  point  sized  particles  with  particles  the  size  of  the 
cell.  No  detail  is  available  at  length  scales  smaller  than  the  cell  size.  Furthermore,  Xd  specifies  the  smallest 
distance  at  which  quasi-neutrality  can  be  assumed.  Motion  of  the  electrons  is  influenced  by  electric  fields 
which  arise  due  to  local  charge  non- neutralities.  A  simulation  cell  which  is  several  orders  of  magnitude  larger 
than  the  Debye  length  is  not  capable  of  resolving  these  charge  variations,  thus  proper  mixing  of  electrons 
with  the  ions  is  not  achieved. 

The  Neumann  (dcjy/dh  =  0)  boundary  condition  was  specified  for  the  potential  solver  on  all  external 
faces.  A  reflective  boundary  condition  was  applied  for  particles  along  the  planes  of  symmetry.  Particle 
motion  through  the  remaining  domain  boundaries  was  controlled  by  an  energy  boundary  condition.  As  was 
described  in  a  greater  detail  in  Ref.  [4],  conservation  of  energy  dictates  that,  in  absence  of  other  potential 
drops  and  energy  sources,  a  particle  traveling  from  rest  through  a  potential  hump  must  reach  a  velocity 
inflection  point.  Since  kinetic  energy  must  remain  non- negative,  the  particle  motion  will  reflect,  and  the 
particle  will  be  trapped  in  a  potential  hill.  Due  to  a  limited  domain  span,  the  inflection  point  may  be 
located  outside  of  the  domain  boundary.  This  case  then  leads  to  a  removal  of  electrons  which  should  have 
been  retained  by  the  simulation.  The  energy  boundary  attempts  to  retain  these  electrons  by  reflecting  any 
particles  with  kinetic  energy  insufficient  to  escape  the  potential  drop.  The  required  potential  energy  is 
calculated  from  the  difference  between  the  beam  core  potential  and  (j)^. 
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C.  Cathode  Model 

Jameson,  et.  al.15  measured  approximately  5V  of  potential  drop  through  the  cathode,  leading  to  electron 
exit  velocities  of  approximately  106m/s.  A  cathode  with  a  keeper  exit  diameter  of  1.2cm  was  used  in  this 
work.  From  j  =  neu,  the  electron  density  at  the  cathode  exit  is  ~  6  x  1016m/s.  Electron  temperature  at  the 
exit  is  ^  leV,  yielding  Ad  ~  3  x  10_5m/s,  or  about  l/10th  of  the  cell  size.  Injection  of  electrons  from  the 
cathode  without  resolving  the  Debye  length  at  the  cathode  tip  resulted  in  development  of  a  virtual  cathode.7 
The  only  electrons  that  were  able  to  escape  were  those  born  with  a  strong  radial  velocity  component. 

In  this  work,  two  cathode  models  are  examined.  The  first  model  uses  a  floating  cathode  with  a  limiting 
value  for  charge  density  near  the  keeper  exit.  This  formulation  is  identical  to  work  done  in  Ref.  [4].  By 
floating  the  potential,  the  strong  axial  potential  gradient  is  eliminated,  and  the  electrons  can  leave  the 
cathode.  Charge  density  around  the  cathode  was  limited  to  prevent  development  of  strong  self-induced 
fields.  The  electrons  were  also  injected  with  small  velocities,  so  they  could  immediately  start  following  the 
electric  field. 


(a)  cathode  subdomain  (b)  number  density 


Figure  3.  The  yellow  shaded  region  indicates  the  subdomain  used  in  cathode  modeling.  Boundary  of  the  mesh 
used  in  neutralization  modeling  and  the  thruster  are  shown  for  scale.  Electrons  were  sampled  from  the  entire 
subdomain  in  subsequent  neutralization  modeling.  Right  plot  (b)  shows  the  electron  density  and  velocity 
stream  lines  at  steady  state. 


While  this  approach  resulted  in  a  free  extraction  of  electrons,  it  also  contributed  to  a  non-physical  plume 
potential  profile  and  electron  temperature.  In  this  work,  a  multi-domain  cathode  model  was  used  to  obtain 
initial  distribution  of  the  neutralizing  electrons.  A  secondary  simulation  was  performed  on  a  sub-domain 
region  surrounding  the  cathode.  Cell  size  was  decreased  to  5  x  10-5m,  and  the  mesh  consisted  of  30  x  40  x  60 
nodes.  Potential  drop  of  10 V  was  applied  between  the  cathode,  and  the  anode,  which  was  represented  by  the 
body  of  the  ion  thruster.  The  subdomain  is  shown  in  figure  3(a).  Boundary  of  the  primary  mesh  is  included 
for  scale. 

Figure  3(b)  shows  the  number  density  of  the  electrons  in  the  sub-domain  after  reaching  a  steady- state. 
Streamlines  of  electron  velocities  are  also  shown.  The  electrons  are  seen  to  expand  uniformly.  Geometry 
of  the  simulation  results  in  a  strong  axial  electric  field  along  the  cathode  centerline.  Electrons  born  in  this 
region  are  accelerated  out  of  the  domain.  The  remaining  electrons  are  slowed  down  by  the  potential  gradient, 
and  are  attracted  back  to  the  anode.  The  velocity  distribution  at  steady-state  is  shown  in  figure  4.  The 
electrons  retain  their  initial  injection  Maxwellian  distribution. 

This  sub-domain  acted  as  a  volumetric  source  during  the  main  neutralization  modeling.  Positions  and 
velocities  of  100,000  randomly  chosen  electrons  were  sampled  to  a  data  file.  A  new  particle  source  has  been 
written  to  return  a  randomly  chosen  entry  from  the  list.  The  potential  drop  between  the  cathode  and  the 
anode  was  fixed  at  10V. 
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Figure  4.  Velocity  distribution  of  the  electrons  in  the  cathode  subdomain. 


IV.  Results 


A.  Reference  Case 


(a)  potential 


(b)  electron  temperature 


(c)  charge  density 


Figure  5.  Potential,  electron  temperature,  and  charge  density  after  3  x  10  7  seconds,  for  a  single  thruster. 

Plots  in  figure  5  show  the  plasma  parameters  on  the  plane  of  symmetry  obtained  using  the  pre-neutralized 
reference  configuration  (RF).  A  distinct  beam  profile  is  seen  in  the  plot  of  the  potential.  Potential  reaches 
about  4.7V  in  the  core  near  the  thruster  exit,  and  is  seen  to  decrease  with  beam  divergence.  This  value 
closely  agrees  with  experimental  measurements.16  Electron  temperature,  computed  by  assuming  Maxwellian 
distribution,  reaches  a  similar  value,  and  is  also  seen  to  decrease  uniformly  with  density.  Figure  5(c)  shows 
the  charge  density,  p  =  q(rii  —  ne).  A  good  neutrality  ratio  is  indicated  by  light-blue  shading.  The  charge 
density  in  the  beam  is  seen  to  be  slightly  positive,  which  leads  to  the  positive  potential  in  the  beam.  An 
electron  sheath  is  seen  to  surround  the  beam.  This  sheath  is  responsible  for  the  containment  of  the  electrons. 

The  electric  field  components,  E  =  —  V0,  are  shown  in  figures  6(a)  and  6(b).  Both  the  radial  and  the 
axial  components  are  approximately  zero  in  the  bulk  of  the  beam.  Hence,  the  acceleration  of  the  electrons 
is  expected  to  be  limited  to  the  regions  near  the  edge  of  the  beam,  with  electrons  moving  at  constant 
velocities  inside  the  beam  core.  The  motion  of  the  electrons  is  highly  random  (fig.  6(c)),  even  though  they 
were  originally  injected  in  the  axial  direction,  using  a  Maxwellian  source  with  Te  =  leV.  Due  to  their  high 
mobility,  the  electrons  seem  to  have  only  a  weak  memory  of  their  injection  distributions. 

Maxwellian  temperature  obtained  from  the  simulation  is  compared  to  the  polytropic  relationship 

/  n  \  fr-1) 

T=T“k)  (8) 
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(a)  electric  field,  x  component 


(b)  electric  field,  z  component 


(c)  electron  velocity  vectors 


Figure  6.  Electric  field  and  electron  velocity  vectors  for  the  reference  case.  Electrons  were  injected  from  the 
optics  using  Maxwellian  distribution  with  Te  =  leV. 


Figure  7.  Comparison  of  numerical  temperature  to  the  polytropic  model,  and  comparison  of  simulation  electron 
density  to  prediction  using  Boltzmann  model. 
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for  three  values  of  7  in  figure  7(a).  Reference  temperature  and  density  were  chosen  to  correspond  to  the  values 
in  the  beam  core,  4.2eV  and  2.5  x  1015  m-3,  respectively.  Neither  of  the  three  chosen  gamma  values  was  able 
to  produce  an  exact  match,  however,  the  temperature  seems  to  roughly  follow  the  polytropic  relationship 
with  7  ~  1.4. 

Numerical  electron  density  was  also  compared  against  the  Boltzmann  relationship,  which  states  that  a 
direct  relationship  exists  between  plasma  potential  and  plasma  density, 

f  (/>-(/>  o\ 
ne  =  n°  6XP  ( ~~kTo~ ) 

Again,  plasma  properties  in  the  beam  core  were  used  for  the  reference  values.  Reference  potential  was  set 
to  4.7V.  The  relationship  was  computed  using  both  constant  reference  temperature  (4.2eV),  and  polytropic 
temperature  with  7  =  1.4.  Generally,  the  agreement  is  not  very  good,  as  figure  7(b)  shows.  Best  agreement  is 
achieved  near  the  core,  which  is  expected,  since  this  location  corresponds  to  the  point  at  which  the  reference 
values  were  sampled.  The  simulation  electron  density  drops  off  faster  than  predicted  by  the  model.  The 
disagreement  is  reduced  by  using  the  polytropic  temperature,  however,  a  significant  discrepancy  still  remains. 
A  better  agreement  could  be  achieved  by  adjusting  the  reference  parameters;  this  approach  however  requires 
prior  knowledge  of  the  solution. 


(9) 


B.  Single  Thruster 

1.  Floating  Cathode 

The  simulation  domain  used  in  single  thruster  neutralization  modeling  can  be  seen  in  figure  3(a).  Due  to 
symmetry,  only  a  half  domain  was  simulated.  Reflective  particle  boundary  condition  was  used  along  the 
symmetric  face.  The  grid  dimensions  were  50  x  100  x  90,  and  the  mesh  used  a  uniform  cell  spacing  of 
2  x  10_4m. 
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(a)  potential 


(b)  electron  temperature 


(c)  charge  density 


Figure  8.  Potential,  electron  temperature,  and  charge  density  after  3  x  10  7  seconds,  for  a  single  thruster. 

Figure  strip  8  shows  potential,  temperature  and  charge  density  obtained  by  floating  the  potential  on  the 
cathode.  Results  for  this  configuration  (NSF)  were  expected  to  agree  with  the  reference  case  RF,  but  a  quick 
comparison  with  figure  8  indicates  a  high  degree  of  divergence.  The  beam  shape  is  no  longer  well  resolved. 
Furthermore,  the  beam  potential  has  increased  to  27V  and  the  maximum  electron  temperature  has  increased 
to  34eV. 

The  high  values  of  beam  potential  and  temperature  point  to  a  non-neutral  beam.  However,  over¬ 
saturation  of  the  beam  with  electrons  by  increasing  the  cathode  current  led  to  turning  of  the  beam  towards 
the  cathode,  but  the  potential  drop  across  the  beam  did  not  change  significantly.  The  beam  temperature 
also  remained  high. 
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(a)  potential 


(b)  electron  temperature 


Figure  9.  Potential,  electron  temperature,  and  charge  density  after  3  x  10  7  seconds,  for  a  single  thruster. 


2.  Multi-domain  Cathode  model 

Therefore,  lack  of  electrons  did  not  seem  to  be  the  main  factor  contributing  to  the  non-physical  results. 
In  this  work,  the  influence  of  the  cathode  model  on  the  results  was  investigated  by  replacing  the  floating 
potential  model  with  the  multi-domain  formulation.  Figure  9  shows  plasma  properties  obtained  with  the 
new  model.  Important  to  note  is  the  region  near  the  cathode.  The  charge  density  plot  shows  a  clear  turning 
of  the  electron  plume  towards  the  beam.  This  behavior  was  not  achieved  in  the  floating  cathode  case.  The 
electron  sheath  surrounding  the  beam  is  also  better  resolved,  and  the  overall  range  of  p  has  decreased  which 
indicates  a  better  neutrality  ratio.  However,  no  significant  improvement  in  results  is  seen.  Although  the 
beam  potential  has  decreased  slightly  to  about  25V,  the  core  temperature  remains  high  and  non-polytropic. 

C.  Analysis 


(a)  RF 


(b)  NSF 


(c)  NSM 


Figure  10.  Electron  density  contour  for  the  three  cases.  A  clear  distinction  is  seen  between  the  reference  case 
and  the  two  cases  in  which  electrons  were  injected  from  the  cathode. 

The  new  cathode  model  failed  to  resolve  the  fundamental  problem  leading  to  the  non-physical  results. 
Figure  10  shows  the  electron  densities  for  the  three  cases.  A  clear  difference  is  seen  between  the  reference  case, 
and  the  two  cathode  cases.  The  electron  density  in  the  reference  case  follows  the  Boltzmann  relationship, 
with  highest  density  coinciding  with  the  beam  core.  The  electrons  instead  seem  to  concentrate  along  the 
edges  of  the  beam  in  the  two  cathode  runs. 
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Existence  of  an  almost  uniform  high-temperature  region  along  beam  axis  in  both  NS  cases  indicates 
mixing  of  electron  streams  with  opposing  directions.  Thus,  the  collective  dynamics  of  the  electrons  seem  to 
be  driven  by  oscillations  around  the  beam  core.  In  other  words,  using  a  1-D  approximation,  electrons  injected 
at  the  cathode  fall  into  the  potential  well  created  by  the  beam.  The  velocity  of  the  electrons  increases  until 
they  pass  the  beam  centerline.  The  velocity  then  begins  to  decrease  as  they  travel  up  the  well.  The  electrons 
come  to  a  stop  at  a  point  where  all  kinetic  energy  has  been  exhausted.  Assuming  initial  injection  at  Om/s, 
the  potential  drop  on  both  sides  of  the  well  will  be  equal.  The  increased  electron  densities  seen  along  the 
beam  edge  in  cases  NSF  and  NSM  are  a  direct  consequence  of  electrons  coming  to  a  stop  and  turning  around. 


(a)  electrons 


Figure  11.  Electron  and  ion  velocity  distribution  at  the  end  of  simulation. 


The  potential  in  the  beam  adjusts  as  a  direct  response  to  the  oscillations  of  the  electrons.  Retention  of 
the  electrons  in  the  beam  requires  that  the  PEbearn  >  KEe^max.  Using  PE  =  KE  and  total  potential  drops 
between  the  cathode  and  the  beam  of  4.7V  and  25V,  the  mean  velocity  of  the  electrons  oscillating  around 
the  beam  core  is  expected  to  be  1.3  x  106  and  2.9  x  106  m/s  for  cases  RF  and  NSx,  respectively.  This  is 
confirmed  by  the  tail  of  electron  velocity  histograms  shown  in  figure  11(a).  The  red  line  (NSF)  shows  the 
distribution  for  the  original  floating  cathode  model.  Comparison  with  the  reference  curve  shows  an  increase 
in  mean  velocity.  Even  more  important  is  the  widening  of  the  curve,  which  is  indicative  of  a  temperature 
increase.  The  shape  of  the  curve  remains  close  to  Maxwellian.  The  green  line  indicates  the  electron  velocity 
obtained  with  the  new  model.  A  small  increase  in  number  of  fast  moving  particles  is  seen,  but  overall,  the 
distribution  remains  comparable  to  NSF.  A  greater  difference  is  seen  in  the  ion  velocities.  Decrease  in  the 
size  of  the  high  potential  region  in  the  NSM  case  results  in  fewer  electrons  being  slowed  down  by  the  strong 
potential  gradient.  The  mean  beam  velocity  for  case  NSM  is  closer  to  the  reference  configuration. 

In  order  to  obtain  the  Boltzmann  relationship,  the  amplitude  of  the  electron  oscillations  must  decrease 
with  time.  This  is  analogous  to  the  classical  example  of  stability,  in  which  a  ball  has  been  placed  into  a 
spherical  cup.  Displacement  of  the  ball  from  the  rest  position  at  the  bottom  will  result  in  simple  oscillations 
about  the  bottom  of  the  cup.  In  the  absence  of  dissipative  forces,  the  amplitude  of  oscillations  will  not 
change  with  time.  However,  in  a  realistic  configuration,  dissipation  of  energy  due  to  non- conservative  forces 
will  result  in  a  gradual  decrease  of  the  amplitude.  After  some  time,  the  ball  will  come  back  to  rest. 

The  simulation  presented  here  does  not  contain  any  such  dissipative  force.  Instead,  the  electrons  keep 
oscillating  around  the  potential  drop  of  the  beam.  Hence,  the  electron  density  is  not  able  to  reach  a 
Boltzmann-like  relationship.  Electrons  in  the  reference  case  are  not  strongly  influenced  by  this  simplification, 
since  they  are  born  at  the  bottom  of  the  potential  well.  Furthermore,  their  initial  velocity  is  coincident  with 
the  direction  of  the  beam  ions.  However,  the  electrons  born  at  the  cathode  originate  at  the  top  of  the 
potential  well  and  flow  into  the  beam  radially. 

Exchange  of  kinetic  energy  with  massive  particles  would  result  in  a  large  decrease  in  electron  velocities, 
while  only  slightly  influencing  the  motion  of  the  ions.  However,  as  was  mentioned  previously,  collisions  were 
ignored  in  this  case,  due  to  large  mean-free  paths,  and  hence  low  collision  frequencies.  More  probable  is  the 
transfer  of  kinetic  energy  from  the  electrons  to  plasma  waves.  Figure  12  shows  the  total  field  energy  versus 
simulation  time.  Clear  oscillations  develop  after  t  =  1.5  x  lO-7  seconds.  The  meaning  of  these  oscillations 
is  still  not  clear,  however  a  decay  is  seen  around  t  =  2.4  x  10-7.  It  is  possible  that  the  numerical  setup  of 
the  problem  is  preventing  development  or  growth  of  energy  dissipating  waves. 
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Figure  12.  Total  simulation  field  energy  versus  simulation  time.  A  high-frequency  wave  is  seen  to  develop 
after  t  =  1.5  X  10— 7  seconds. 

V.  Conclusion 

A  new  simulation  model  to  study  ion  beam  neutralization  is  being  developed.  This  model  uses  a  fully- 
kinetic  formulation  in  which  both  electrons  and  ions  are  tracked  as  macro-particles.  This  formulation  avoids 
problems  associated  with  fluid  modeling  of  the  electrons,  but  introduces  numerical  difficulties.  Most  im¬ 
portant  is  the  necessity  to  resolve  the  local  Debye  length,  otherwise  the  electrons  fail  to  mix  with  the  ions. 
Computationally  excessive  number  of  nodes  would  be  required  to  resolve  the  Debye  length  on  a  full-scale 
geometry.  Instead,  a  dimensional  scaling  approach  was  used,  with  thruster  current  adjusted  such  that  plasma 
density  at  the  thruster  exit  remained  unchanged.  This  approach  allowed  for  examining  the  neutralization 
process.  Out  flux  of  electrons  at  boundaries  was  prevented  by  reflecting  all  electrons  with  kinetic  energies 
insufficient  to  escape  the  potential  drop  of  the  beam. 

This  simulation  approach  was  used  to  model  neutralization  of  the  NASA  NEXT  ions  thruster.  A  reference 
case  was  setup  by  injecting  both  electrons  and  ions  from  the  optics.  The  potential  solution  showed  a  clear 
beam  profile,  with  maximum  potential  of  4.7V.  The  electron  temperature  reached  about  5eV  in  the  core, 
and  decreased  polytropically  with  density.  These  results  agree  well  with  experimental  data.  The  electron 
density  was  also  compared  to  the  Boltzmann  model,  but  the  two  curves  diverged  for  the  chosen  coefficients. 

Simulation  of  a  single  thruster  neutralized  by  an  external  cathode  was  also  studied.  The  potential  around 
the  cathode  could  not  be  resolved  correctly  using  the  primary  mesh,  due  to  the  cathode’s  small  size  and 
a  high  electron  density  near  the  tip.  Instead,  two  approaches  were  investigated.  In  the  first  model,  the 
cathode  potential  was  allowed  to  float  and  charge  density  at  the  tip  was  fixed.  The  second  model  used  a 
multi-domain  formulation.  Simulation  of  the  near  cathode  region  was  performed  first  using  a  fine  mesh, 
followed  by  sampling  of  random  electrons  to  a  data  file.  This  distribution  list  was  then  used  by  primary 
simulation  to  introduce  electrons  from  the  volume  described  by  the  fine  mesh. 

Plasma  properties  in  either  cathode  run  did  not  agree  with  the  reference  case.  Beam  potential  increased  to 
about  25V  and  the  temperature  distribution  no  longer  followed  the  polytropic  relationship.  Closer  inspection 
of  simulation  results  indicates  that  the  electrons  oscillate  around  beam.  The  amplitude  of  oscillations  should 
decrease  with  time,  however,  this  does  not  seem  to  be  the  case  in  the  current  simulation.  Investigation  into 
the  primary  mechanism  and  the  actual  modeling  of  the  kinetic  energy  sink  remain  as  future  work. 
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